A solution of (\ref{eqn:waveeqnsept}) is assumed to have the form,
\begin{align}
		T(t) = e^{{st}}\label{eqn:solnT}
\end{align}
Whose first and second derivatives are the following,
\begin{align}
		T'(t)&={s}e^{{st}}\label{eqn:dsolnT}\\
		T''(t)&={s}^2e^{{st}}\label{eqn:ddsolnT}
\end{align}
Substituting (\ref{eqn:solnT})--(\ref{eqn:ddsolnT}) into (\ref{eqn:waveeqnsept}) yields,
\begin{align}
\mu\varepsilon{s}^2e^{{st}}+\left(\mu\sigma_e+\varepsilon\sigma_m\right)se^{{st}}+(\sigma_m\sigma_e-\gamma^2)e^{{st}}&=0\label{eqn:waveeqnseptsoln1}
\end{align}
which simplifies to the following,
\begin{align}
\left[\mu\varepsilon{s}^2+\left(\mu\sigma_e+\varepsilon\sigma_m\right)s+(\sigma_m\sigma_e-\gamma^2)\right]T(t)&=0\label{eqn:waveeqnseptsoln2}
\end{align}
Since $T(t)$ is never zero for real values of $s$, it is apparent the only way to satisfy the differential equation is to satisfy the following quadratic equation,
\begin{align}
\mu\varepsilon{s}^2+\left(\mu\sigma_e+\varepsilon\sigma_m\right)s+(\sigma_m\sigma_e-\gamma^2)&=0\label{eqn:waveeqnseptsoln3}
\end{align}
This equation is called the \emph{characteristic equation (Eigen equation)}, whose solutions are called the \emph{characteristic roots (Eigen values)}.  In general the quadratic will have two solutions $s_1$ and $s_2$.  
\begin{align}
s_{1,2}&=-\frac{1}{2}\left(\frac{\sigma_m}{\mu}+\frac{\sigma_e}{\varepsilon}\right)\pm\sqrt{\left[\frac{1}{2}\left(\frac{\sigma_m}{\mu}+\frac{\sigma_e}{\varepsilon}\right)\right]^2+\frac{\gamma^2}{\mu\varepsilon}-\frac{\sigma_m}{\mu}\frac{\sigma_e}{\varepsilon}}
\end{align}
There are three possible cases of solutions for the characteristic equation, distinct real roots, real equal roots, and a complex conjugate pair of roots if $T(t)$ is assumed to be a real valued function.

\subsubsection{Distinct Real Roots}
Under the assumption that that the characteristic equation has two unequal real roots $s_1$ and $s_2$, there are two solutions,
\begin{align}
	T_1(t) &= c_{1}e^{s_{1}t}\\
	T_2(t) &= c_{2}e^{s_{2}t}
\end{align}
These functions are linearly independent, therefore they form a fundamental set. This can be verified that the \emph{Wronskian} is not equal to zero.  The general solution is obtained by superposition.
\begin{align}
	T(t) = c_{1}e^{s_{1}t}+c_{2}e^{s_{2}t}\label{eqn:DistinctRealRoots}
\end{align}
\subsubsection{Repeated Real Roots}
When $s_1=s_2$ the solution to the differential equation is the following,
\begin{align}
	T(t) = c_{1}e^{s_{1}t}+c_{2}te^{s_{1}t}\label{eqn:RepeatedRealRoots}
\end{align}
As long as $s_1<0$ the signal will approach zero as $t$ approaches infinity.
\subsubsection{Complex Conjugate Roots}
If the roots are complex conjugates then,
\begin{align}
 s_1 &= \sigma+j\omega\\
 s_2 &= \sigma-j\omega
\end{align}
where $\sigma>0$ and $\omega>0$. Formally, there is no difference between this case and the case of distinct real roots, therefore,
\begin{align}
	T(t) = c_{1}e^{(\sigma+j\omega)t}+c_{2}e^{(\sigma-j\omega)t}\label{eqn:solnDistinctComplexRoots}
\end{align}
By using Euler's formula,
\begin{align}
	e^{j\theta}=cos\theta+j\sin\theta
\end{align}
equation (\ref{eqn:solnDistinctComplexRoots}) can be rewritten as,
\begin{align}
	T(t) &= c'_{1}e^{\sigma{}t}\cos(\omega{}t)+c'_{2}e^{\sigma{}t}\sin(\omega{}t)\label{eqn:ComplexConjugateRoots}\\
	&= e^{\sigma{}t}\left(c'_{1}\cos(\omega{}t)+c'_{2}\sin(\omega{}t)\right)\label{eqn:ComplexConjugateRoots2}
\end{align}
where $c'_1 = c_1+c_2$ and $c'_2 = j(c_1-c_2)$.

\subsubsection{Time Harmonic Solutions}
To obtain time harmonic solutions the condition of $\sigma=0$ in (\ref{eqn:solnDistinctComplexRoots}) must be enforced.  Either solution in (\ref{eqn:solnDistinctComplexRoots}) is a valid solution and since they are assumed causal, only one of them is needed.  Therefore it is necessary to choose one solution as a standard and use it consistently.  The solution chosen will be,
\begin{align}
	T(t) &= e^{j\omega t}\label{eqn:steadystate}
\end{align}
where its derivative is given as,
\begin{align}
	T'(t) &= j\omega e^{j\omega t}\label{eqn:dsteadystate}
\end{align}
The result in (\ref{eqn:steadystate}) was obtained by letting,
\begin{align}
	s=j\omega\label{eqn:sstandard}
\end{align}
Substituting (\ref{eqn:sstandard}) into (\ref{eqn:waveeqnseptsoln3}) yields,
\begin{align}
\gamma^2=-k^2=-\omega^2\mu\varepsilon+\sigma_{m}\sigma_{e}+j\omega(\mu\sigma_{e}+\varepsilon\sigma_{m})\label{eqn:gk0}
\end{align}
This can be rewritten as,
\begin{align}
\begin{split}
\gamma^2=-k^2&=\omega^2\mu\varepsilon\left(-1+\sigma_m'\sigma_e'+j(\sigma_m'+\sigma_e')\right)\\
&=-\omega^2\mu\varepsilon(1-j\sigma_m')(1-j\sigma_e')\\
&=-\omega^2\mu'\varepsilon'
\end{split}\label{eqn:gk1}
\end{align}
where,
\begin{align}
\sigma_m' &= \frac{\sigma_m}{\omega\mu}\label{eqn:sigma_m_primed}\\
\sigma_e' &= \frac{\sigma_e}{\omega\varepsilon}\label{eqn:sigma_e_primed}
\end{align}
and where,
\begin{align}
\mu' &= \mu(1-j\sigma_m')\label{eqn:mu_complex}\\
\varepsilon' &=  \varepsilon(1-j\sigma_e')\label{eqn:eps_complex}
\end{align}
The variable $\gamma$ is in general complex and the real and imaginary parts can be explicitly solved for.  Let $\gamma=\alpha+j\beta$ and substitute into (\ref{eqn:gk1}) to get,
\begin{align}
(\alpha+j\beta)^2=(\alpha^2-\beta^2)+j2\alpha\beta=\omega^2\mu\varepsilon\left(-1+\sigma_m'\sigma_e'+j(\sigma_m'+\sigma_e')\right)
\end{align}
Equating real and imaginary parts yields,
\begin{align}
\alpha^2-\beta^2&=\omega^2\mu\varepsilon^2(\sigma_m'\sigma_e'-1)=A\\
2\alpha\beta&=\omega^2\mu\varepsilon(\sigma_m'+\sigma_e')=B
\end{align}
Solving these equations, they can be manipulated to the following form,
\begin{align}
\alpha^4-A\alpha^2-\left(\frac{B}{2}\right)^2&=0\\
\beta^4+A\beta^2-\left(\frac{B}{2}\right)^2&=0
\end{align}
The solution to each of these have in general four roots, however only the roots that yield positive values for $\alpha$ and $\beta$ are chosen. Solving for $\alpha$ and $\beta$ gives,
\begin{align}
\alpha&=\frac{1}{\sqrt{2}}\sqrt{\sqrt{A^2+B^2}+A}\\
\beta&=\frac{1}{\sqrt{2}}\sqrt{\sqrt{A^2+B^2}-A}
\end{align}
Finally, substituting in the constitutive parameters yields,
\begin{align}
\alpha&=\frac{\omega\sqrt{\mu\varepsilon}}{\sqrt{2}}\sqrt{\sqrt{1+\sigma_m'^2+\sigma_e'^2+(\sigma_m'\sigma_e')^2}+\sigma_m'\sigma_e'-1}\label{eqn:alpha}\\
\beta&=\frac{\omega\sqrt{\mu\varepsilon}}{\sqrt{2}}\sqrt{\sqrt{1+\sigma_m'^2+\sigma_e'^2+(\sigma_m'\sigma_e')^2}-\sigma_m'\sigma_e'+1}\label{eqn:beta}
\end{align}
When $\sigma_m=0$ then (\ref{eqn:alpha}) and (\ref{eqn:beta}) reduce to the following.
\begin{align}
\alpha&=\frac{\omega\sqrt{\mu\varepsilon}}{\sqrt{2}}\sqrt{\sqrt{1+\sigma_e'^2}-1}\label{eqn:alpha1}\\
\beta&=\frac{\omega\sqrt{\mu\varepsilon}}{\sqrt{2}}\sqrt{\sqrt{1+\sigma_e'^2}+1}\label{eqn:beta1}
\end{align}
When $\sigma_m=\sigma_e=0$ then (\ref{eqn:alpha}) and (\ref{eqn:beta}) reduce to the following.
\begin{align}
\alpha&=0\label{eqn:alpha2}\\
\beta&=\omega\sqrt{\mu\varepsilon}\label{eqn:beta2}
\end{align}

\subsubsection {Higher Order Media}
Materials that have constitutive parameters that are a function of frequency are characterized as \emph{non-linear}.  One possible way to model such materials is to add higher order temporal derivatives into Maxwell's equations.  Just as the fictitious magnetic loss $\sigma_m$ was added to Maxwell's equations to balance the equations, higher order derivatives with constant coefficients will be added to model the non-linear properties.  Consider the following set of modified Maxwell's equations,
\begin{align}
    \nabla\cdot\vec{\mathcal{E}}&=0&&\text{Gauss's Electric Field Law}\label{eqn:geflmod}\\
    \nabla\cdot\vec{\mathcal{H}}&=0&&\text{Gauss's Magnetic Field Law}\label{eqn:gmflmod}\\
    \nabla\times\;\vec{\mathcal{E}}&=-\sum_{m=0}^{M}a_m\frac{\partial^m\vec{\mathcal{H}}}{\partial{t}^m}&&\text{Faraday's Law}&&\label{eqn:flmod}\\
    \nabla\times\,\vec{\mathcal{H}}&=\sum_{n=0}^{N}b_n\frac{\partial^n\vec{\mathcal{E}}}{\partial{t}^n}&&\text{Ampere's Law}\label{eqn:almod}
\end{align}
Taking the curl of both sides of (\ref{eqn:flmod}) yields,
\begin{align}
\nabla\times\nabla\times\vec{\mathcal{E}}&=-\sum_{m=0}^{M}a_m\frac{\partial^m}{\partial{t}^m}\nabla\times\vec{\mathcal{H}}\label{eqn:curlflsfmod}
\end{align}
Using the vector identity $\nabla\times\nabla\times\vec{A} = -\nabla^2\vec{A}+\nabla(\nabla\cdot\vec{A})$ and substituting (\ref{eqn:almod}) into (\ref{eqn:curlflsfmod}) gives,
\begin{align}
-\nabla^2\vec{\mathcal{E}}+\nabla(\nabla\cdot\vec{\mathcal{E}})&=-\sum_{m=0}^{M}a_m\frac{\partial^m}{\partial{t}^m}\sum_{n=0}^{N}b_n\frac{\partial^n\vec{\mathcal{E}}}{\partial{t}^n}\label{eqn:waveeqn1mod}
\end{align}
Rearranging and using the fact that the divergence is equal to zero, (\ref{eqn:geflmod}) gives the following wave equation,
\begin{align}
\nabla^2\vec{\mathcal{E}}&=\sum_{m=0}^{M}\sum_{n=0}^{N}a_mb_n\frac{\partial^{m+n}\vec{\mathcal{E}}}{\partial{t}^{m+n}}\label{eqn:waveeqn2mod}
\end{align}
By following a similar procedure a wave equation for the magnetic field intensity can also be derived.  Letting $M=N$, this can be written in a more compact format as,
\begin{align}
\nabla^2\vec{\mathcal{E}}&=\sum_{p=0}^{2N}c_p\frac{\partial^{p}\vec{\mathcal{E}}}{\partial{t}^{p}}\label{eqn:waveeqn3mod}
\end{align}
where,
\begin{align}
C_p=
\begin{cases}
\sum_{m=0}^{p}a_{m}b_{p-m}\;\;&\text{if}\;\;  0\le p \le N\\
\sum_{m=p-N}^{N}a_{m}b_{p-m}\;\;&\text{if}\;\;  N+1\le p \le 2N
\end{cases}
\end{align}
